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Abstract 

We consider the effect of correlated noise on the stability of synchronisation 
of oscillators on a general network. By examining the Kuramoto model in 
the neighborhood of a global phase synchronised fixed point the impact of 
the noise is seen in the time-dependent probability density. If the support of 
the distribution evolves outside the basin of attraction of a fixed point with 
finite probability the system is deemed to be unstable. By exactly solving the 
Fokker-Planck equation we find that quite different instabilities follow. For 
uncorrelated noise the instability is exponentially suppressed: for small dif- 
fusion constant, there is a vanishingly small probability that the system can 
drift out of phase synchronicity. With correlated noise applied to oscillator 
frequencies and the network coupling, the peak of the probability distribu- 
tion itself drifts outside the basin of the fixed point and suppression becomes 
power-law. Correlated noise can therefore strongly de-stabilise phase syn- 
chronicity. The significance of result for general networks is discussed. 
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1. Introduction 

The spontaneous generation of ordered or patterned behaviour is a sig- 
nificant property of nonlinear dynamical systems of many coupled heteroge- 
neous entities. This behaviour is argued to be relevant to a vast array of real 
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world systems: physical, chemical, biological, ecological and social. Where 
the dynamics of the component level entities have cyclic characteristics, the 
model of coupled phase oscillators on a network is a compact mathematical 
representation offering much insight. The Kuramoto model [ij] is one espe- 
cially simple example of this, exhibiting a clear transition from incoherence 
to phase synchrony as a single coupling constant is varied from small to large 
positive values. In this paper we explore how different ways of injecting noise 
- or stochastic fluctuations - from a regime close to the synchronised state 
can generate instability. 

The Kuramoto model [l[ involves one-dimensional phase oscillators de- 
scribed by angles 9i coupled on a network described by an undirected graph 
G of size and is given by the first order differential equation, 

K ^ 

= - ^ XI -Oj) ,l<i<N. (1) 

i=i 

Here Aij is the adjacency matrix encoding the graph structure, Ui are in- 
trinsic frequencies of the oscillators, usually drawn from some statistical dis- 
tribution, and is a global coupling controlling the rapidity and strength 
of mutual adjustment between adjacent oscillators. A good overview of this 
and related models is given by Strogatz At zero coupling the system 
behaves incoherently with each oscillator moving in simple harmonic motion 
according to its intrinsic frequency. Increasing the coupling drives the popu- 
lation into synchrony, with some of the oscillators locking into a single core 
that moves with harmonic motion according to the frequency average over 
the whole oscillator ensemble. Kuramoto jl| originally considered only the 
complete {Aij = 1) infinite (A^ — )■ oo) case, which lends itself to solution via 
a variety of techniques from fluid dynamics. In this context, stability and 
noise have been studied using the continuum limit js], 01 . However, these 
techniques are difficult to generalise to non-complete finite graphs. Numeri- 
cal studies, such as jH, have shed light on the path towards synchronisation 
for different network topologies but it is difficult to draw specific conclusions 
about stability from these. 

A contrasting class of dynamical models of the form 

0i = fi<P^) + Mai^j) - am (2) 
j 

have also been tested for stability in the vicinity of fixed points characterising 
synchrony. Using spectral graph theory, Pecora and Carroll ^ formulated a 
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Master Stability Equation from Eq. which identified stabihty as a property 
of the spectrum of the graph Laplacian (to be defined below) of the underly- 
ing network. As usual in the Lyapunov approach, a fixed point (^* (associated 
here with phase synchrony) is perturbed by fiuctuations 0* + rji and inserted 
into the dynamical equations. By isolating powers of rj up to first order {rj is 
'small': 77^ ~ 0) solutions for both the fixed point and the fiuctuations can be 
obtained. In this respect, the difference between the Kuramoto model Eq.(IT]) 
and Eq.([2]) lies in the interactions of the former depending on a function 
of phase differences and the latter a difference of a function of the phases. 
This seemingly innocuous distinction means that within such first order con- 
siderations the phase synchrony of the Kuramoto model is globally stable 
(one must go to second order in fiuctuations to detect instabilities here j3]) 
while for models of the form Eq.([2]) rich instabilities can be identified more 
straightforwardly. 

However, not all properties of real world systems, such as the 'fine- 
structure' of component entities below the length or time scales of human 
observation, can be systematically encoded in deterministic equations of mo- 
tion. This is where 'noise' or stochastic fiuctuations in equations of motion 
are invoked. The stability in the presence of noise of a key nonlinear sys- 
tem property such as synchrony is therefore important in eventually applying 
such elementary mathematical models to real world systems. Park and Kim 
js], lof have numerically examined noise in the Kuramoto model but for the 
limited case of equal frequencies and a complete network (and also using an 
extra 'pinning term'). While giving deep insights into new types of meta- 
stable states for this case, no conclusions can be drawn for inhomogeneous 
frequency distributions and general graphs (each case of interest has to be 
separately simulated). 

The ability to treat general graphs and frequency distributions is where 
analytic approaches, such as the Lyapunov approach to deterministic stabil- 
ity, come into their own. However going beyond the Lyapunov approach for 
the Kuramoto model has its challenges. Fluctuations for the deterministic 
case correspond to instantaneously shifting the system from a fixed point but 
then allowing the deterministic system behaviour of Eqs. (ITPj) govern the sub- 
sequent evolution. Stochastic fiuctuations, contrastingly, involve sustained 
erratic perturbation away from the deterministic trajectories. Mathemati- 
cally this is not so easily cast into the fiuctuation equations like rj (x rj oi 
deterministic systems. For the Kuramoto model the nonlinearities are non- 
polynomial (a cosine rather than, say, quartic potential energies whose turn- 
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ing points define a boundary across which 'first passage time' for diffusion can 
be computed). A shghtly weaker approach solves the Fokker-Planck equation 
for the time-dependent probabihty distribution for noise fluctuations in the 
vicinity of a fixed point. One identifies then whether the stochastic dynamics 



will keep the system in that vicinity or allow it to escape [10[. Instability of 
the stochastic system can be thereby detected. 

In terms of the Kuramoto model in the vicinity of phase synchrony noise 
be injected by three basic ways: through the coupling, the frequencies, and 
the network topology. We consider combinations of the first two methods 
while treating the underlying network fixed. However, in the neighborhood 
of phase synchrony the system decouples into independent behaviours of 
modes corresponding to eigenvectors of the graph Laplacian just as in the 
approach of Pecora and Carroll [cj to systems like Eq. ([2]) . We therefore apply 
the noise (to coupling and frequencies) independently for each such mode of 
the system. When this noise is introduced correlated strong stochastic in- 
stabilities are generated, in the sense just discussed. We solve, for a variety 
of cases, the stationary probability dsitributions. In some cases we are able 
to go even further and analytically solve for the time-dependent probability 
distributions. In considering the application to general networks we are able 
to deduce that phase synchronisation on networks with many poorly con- 
nected subgraphs (such as sparse random graphs and scale-free graphs) is 
more susceptible to stochastic instability. 

The paper is structured as follows. We formulate the model in the next 
section, and then discuss our solution method for the Fokker-Planck system. 
To streamline the paper we next state the main results for the probability 
distribution functions with minimal derivation (the laborious detail is laid 
out in Appendices at the end of the paper). We then deduce the cases where 
instability applied. Subsequently we consider some special parameter cases 
and plot the probability density functions. The paper concludes with a dis- 
cussion where we draw implications of our results for general graphs, compare 
our results to known cases in the literature and look ahead to future work. 
The detail of our derivations is presented in a series of seven appendices. 
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2. Stochastic equations close to synchrony 



2.1. Laplacian decomposition of deterministic case 

We wish to consider the system Eq.([T]) for a connected graph G close to 



synchrony. To this end it is useful to use the Laplacian for G, 

Lij = diSij ^ij; (3) 
with di is the degree of node i. The spectrum of Laplacian eigenvalues Aj 



is positive semi-definite, as explained in [Appendix A[ The corresponding 



system of orthonormal eigenvectors is given by, 

N N 
j=l i=l 

Some relevant theorems involving the Laplacian spectrum are summarised in 



Appendix A Most immediately useful here is the property that the lowest 



eigenvalue vanishes, Aq = 0. For a connected graph G the first non-zero 



eigenvalue Ai is positive and known as the algebraic connectivity [12|. The 
corresponding eigenvector is called the Fiedler. Graphs with Ai < 1 are 
easily disconnected by removal of a small number of links. The eigenvectors 
z/^'') also have a network interpretation. They represent weighted sub-graphs 
of G: the (not necessarily positive) numerical weight with which a node i is 

(r) 

represented in the Laplacian eigenmode r is given by the component . 

The Laplacian is in fact a Laplace-Beltrami operator on the manifold de- 
scribed by G. To this extent it reflects properties well-known in condensed 
matter systems: the eigenvalues A^ are a 'momentum- squared' and their cor- 
responding eigenvectors act like 'probes' on substructures of the graph. The 
eigenvectors corresponding to low-lying eigenvalues r ~ 1 involve larger scale 
sub-graphs of G while those for higher values r ^ N reflect finer structures 
in the graph. This will be useful in interpreting our results. 

We now expand a phase angle of the Kuramoto model in the Laplacian 
basis of eigenmodes, 

Af-l 

^,,(t) = ^x.(t)z.f\ (4) 

r=0 

Treating Ui as the components to an A^-dimensional vector, cu, with u the 
average over the frequency distribution, we can form the scalar product of 



5 



the frequency vector and the Laplacian eigenvectors 



Turning to the Kuramoto model, and combining K/N into a couphng con- 
stant cr, we consider the system close to a point of global phase synchronisa- 
tion [Oi ^ 6j 

Qi ^ cui - o- AijiQi - Oj) 
j 

= Ui- aJ^LijOj. (5) 



Inserting the expansion Eq.(jl]) into the truncated equation of motion 
Eq-dS]), we can extract equations for the normal eigenmodes r 7^ 0, 

Xr{t) = W*-''-* — aXrXr{t) (6) 



while the zero mode of Eq.([S]) gives Xo{t) = yN{t — t')uj at time t after 
initial time t' . To this order, the modes r are decoupled. The exact solution 
to Eq.® is: 

Xr{t) = x^e''^"'-^*-*') + ^ f 1 - e-'^"'-(*-*')) (7) 

where x'^ = Xr{t') represents the initial configuration. Taking the large t limit 
of this shows that there is a fixed point 

< = V (8) 

whose stability is determined by the exponents —cxXr. These are effectively 
the Lyapunov exponents and they do not change sign due to the semi- 
definiteness of the Laplacian spectrum (see Appendix A ). As discussed by 
one of us elsewhere [7|], this shows that global phase synchrony is stable in 
every direction of trajectory space (which can be labelled by the modes r): 
fluctuations on 9i = 9j \/ i,j are small, 

|x,(t)|<lVt (9) 
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as long as the condition is true at the initial time t = t', and 

|a;;| = | — 1<1 (10) 

SO that the approximations leading to Eq.Q are valid over all times. The 
quantity x* defines the origin of a 'basin of attraction' in trajectory space 
whose curvature in every direction of the origin is concave (since the aXr do 
not change sign). 

This leads to the following simple intuitive picture to which we shall 
often return. The Laplacian eigenmodes with largest values of indi- 
cate whether the system can phase synchronise or not. This ratio brings 
together the three key elements of the dynamical system: the frequencies, 
the network topology (in the eigenvalues and eigenvectors) and the coupling 
constant. There are a number of ways these may combine to satisfy Eq.f llOp . 
Two obvious cases are sufficiently strong coupling and/or sufficiently narrow 
frequency range. Modes r with large also satisfy Eq. lfTOj) . Easily discon- 
nected graphs have < 1 and therefore require larger coupling to globally 
phase synchronise. This is also intuitively sensible. However - in light of our 
comments about eigenvectors above - the A^ for r ~ 1 correspond to coarser 
structures so that there will be more nodes represented in the eigenvector 
giving possibly larger values of uj^'^\ Thus the mode most susceptible to 

(r) 

disruption from synchronicity is that for which is closest to one. We call 
this the critical mode. Because of the interplay of quantities in the ratio in 
X* this mode may not be that corresponding to r = 1. 



We obtain a weak type of criterion for synchronisation (see also 13|). 
a ^ \uj^'^^\/ Xr^r. But because we cannot determine a boundary to the basin 
we therefore cannot, at this order, extract a sharp threshold of coupling 
which distinguishes synchrony from incoherence as distinct phases. To de- 
tect instability for the deterministic case one must go to second order in 
fluctuations [7|. Possible instabilities induced by noise, the stochastic case, 
are our concern in this paper. 

2.2. Stochastic system in Laplacian decomposition 

The famous logistic equation x = x{l — x) for population dynamics can 
be brought to the form Eq.(|6]) after a change of variables. This is our clue 
for examining this system in the presence of stochastic fluctuations. In terms 
of phase synchronisation, this relationship is understood by seeing the nor- 
mal mode fluctuations Xr in terms of population species j7|. The Laplacian 
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normal modes correspond to a variety of 'species'. Extinction of a species 
amounts to approach to the phase synchronised fixed point, namely stability 
of the fixed point; extinction of all species is equivalent to the system be- 
ing at the phase synchronised fixed point. Growth corresponds to departure 
from the phase sychronised fixed point, namely instability of the fixed point. 
Bringing stochastic fluctuations into the picture now enables us to exploit 
a significant body of work on noisy population dynamics. Specifically, in 



14( 1 and references therein, Gofa considered a generalisation of the logistic 



equation by perturbation with two correlated Gaussian White Noise (GWN) 
terms, vi""^ and vi^^; the superscripts a and m stand for additive and mul- 
tiplicative respectively. We perform a similar analysis by perturbing Eq.(l6]) 
for each eigenmode via the following substitutions, 

(r) 

a ^ a+^ri^\t) (11) 
where |7i^\72^''| G IR- The resulting Langevin equations are 

Xr{t) = (u'-'^ + l^2^rl!^\t)) - {aXr + li'^T^r^t)) x,(t) (12) 

We see the meaning now of the terms 'additive' and 'multiplicative' in terms 
of how they manifest noise into the equation of motion: the former introduce 
a time dependent randomness on the intrinsic frequencies, u^^\ and the latter 
term introduces it on the coupling, a. We emphasise that the fluctuations 
Eq.( lTT|) are introduced for each Laplacian eigenmode, leading to a localisation 
of the effective coupling. 

In light of this, we adopt an approach which probes specifically the sen- 
sitivity of the aforementioned critical modes to noise by targeting stochastic 
fluctuations to specific clusters of the network. The clusters are determined 
by the Laplacian eigenmode decomposition. 

The additive and multiplicately GWN terms are themselves given in terms 
of uncorrelated GWN functions vi^^ and vi^^ through 



Ti^\t) = CrTl'^t) + v^T^r(2)(t) , r^'") = r« , g m, (13) 

where the following expectation values encode the absence of correlation 

2). 

(r«(t)) = (r(2)(t)) = (r«(t)rf(0) = o 



between vi^^ and Fr : 
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(rW(t)rSnt')) = (n?Wr(?(t')) = nSr,rAt-t'). 

The quantity Q is the diffusion constant. As Q approaches zero we expect 
to obtain the usual deterministic results. Hence vi and Fr™'' are correlated 
in the following way, 

(r(^)(t)r(™)(t')) = Cr,n5r,rAt -t'),neR (i4) 

We refer to the constant Cr as the correlation parameter. For more detail on 
Langevin equations and GWN the reader is referred to Chap. 3 of [l5|. 

Rather than focusing on Xr in the Langevin equations (or specifically, 
associated 'random variables' X^) we shall address the corresponding proba- 
bility density function associated. Denoted P, it determines the probability 
that the system lies in a state Xr < Xr < Xr + dxr- The Fokker-Planck 
equation governs the evolution of this probability density, given the system 
is initially in state X^. = x'^ at initial time t': 

P{x,t) = P{x,x'\t,t') 
P{x,t') = 6{x-x') 

where x now represents the vector of components Xr- From the Langevin 
equation Eq. f|T2|) we obtain, using 15|], the following Fokker Planck equation, 

= E 9^ {s^'\xr)Pix,t)} - 5^ — {g«(x.) P{x,t)} (15) 



r=l r=l 



where, 

sW(x^) 







^ < 




V) _ 





Cr e [-1, 1] 
1 1 c, gM/[-i,i] 



(16) 

Here s^^^ represents a diffusion matrix and drift vector. In the form 

Eq.f lT^ . the diffusion matrix has been diagonalised. 



3. Solving the Fokker-Planck equations 

3.1. Outline of solution method 

In the following, we sketch out the method of solving the Fokker-Planck 
equations in order to present relatively quickly the final solutions. Detailed 
derivations are given in the appendices. 
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Due to the diagonal form of the diffusion matrix, we decompose P{x, t) 
into a product of — 1 terms, 



N-l 



P{x,t) = Yl Pr{Xr,t) 



r=l 

SO that Eq. flT^ decouples into separate equations for each r 

-PriXr,t) = —{s^^)(x)PriXr,t)}- — {q^'\x)PriXr,t)} (17) 

forr = {l,...,A^ — 1}. Since we have decoupled the modes we shall drop the 
sub (super) scripts r from this point on. Solving Eq. ([T7|) is the main concern 
for the remainder of this paper. 

3.2. Stationary densities 

We first solve for the stationary density of Eq. (fT7|) . denoted Pst{x), by 
setting -^P(x,t) = 0. We obtain the corresponding (Pearson's) differential 
equation, 

-^{.s{x)Pstix)}-q{x) Pst(x) = 0. 
ax 

Being linear its solution is 

P..(,) = _i_exp{/<;xfM}. (18) 

Pst is also referred to as a weight function. Substituting the product 

P{x,t) = Pst{x)g{x,t) (19) 
into Eq.l fTTl) . the ensuing equation for g{x,t) is 

d d"^ d 

— ^(x, t) = s{x)—g{x, t) + q{x)—g{x, t). (20) 
We identify a Sturm-Liouville operator 
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The factor g{x,t) can thus be decomposed in terms of the normal modes of 

n 

nu^ix)=^u^{x),^^0. (22) 

The eigenvalues ^ may be discrete or continuous depending on boundary 
conditions, to be specified exphcitly below. Then the full probability dis- 
tribution can be expressed, as usual, in terms of these eigenfunctions and 
eigenvalues 

p{x,x'\t,t') = p,t{x)J2p^(^'^^''''^MxMi^') 

where the sum will be an integral for the continuous parts of the spectrum, 
and constants should be determined so that the initial condition at t = t' 

P{x,x'\t',t') = 6{x-x') 

is satisfied. 



3.3. Exact solutions to Fokker-Planck equations 



In contrast to [IJI , we can go somewhat further than the stationary parts 
and solve - in some cases - the time dependence. Here we lay out the structure 
of that solution for the cases we have successfully solved. 

Determining the nature of the spectrum of eigenvalues ^ and eigenfunc- 
tions are the boundary conditions in x, x' - oscillatory or non-oscillatory 
- which in turn relate to the order of the coefficient polynomial s{x). From 
Eq. (fT6!) . s{x) is either a constant {'^[^^ = 0) or a quadratic polynomial 
{-f[^^ ^ 0). Three cases of boundary conditions are relevant: 

• Case I: Non-oscillatory at boundaries, s{x) constant: the spectrum 
is entirely discrete and the eigenfunctions u are Hermite polynomials. 
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Case II: Oscillatory at one boundary, s{x) quadratic: the spectrum is 
mixed, with the discrete eigenvalues associated with Bessel polynomials 



17| as eigenfunctions u and eigenfunctions for the continuous spectrum 



are the Whittaker functions. 
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• Case III: Oscillatory at both boundaries, s{x) quadratic: the spectrum 
is again mixed. For this case, however, we have not been able to satis- 
factorily solve the full time-dependence (although we know the eigen- 
functions u for the discrete spectrum are the Romanovski polynomials 
fisf). We shall only give the stationary density for this case. 

So that the reader can readily see the overall form of the solutions, we 
present them here without detailed derivation. 

3.3.1. Case I 

In terms of a variable z = a/ aX/^l{x — u/aX), the probability density 
function for this case is given in terms of Hermite polynomials Hn{z): 



P^'\z,z%t') = ^^e-^^$: J-^e--"(-')/7„(z)if„(.') (23) 
The prefactor originates from the stationary part of the distribution. 



After applying the Mehler formula for Hermite polynomials (see Appendix F ) 
the final form of the probability density becomes. 



P(') {x,x'\t,t') 



a 



A r a\ 



7r7|fi(l - e-2-A(i-i'))' "^""P 1^ ^2^(1 _ e-2-A(i-t')) 

X (x - x'e-'^^(*-*') - ^ (l - e-'^^(*-*')))'| . (25) 

We recognise this as a Gaussian smearing around the deterministic solu- 
tion, Eq.(l7]). In the appendix we show that for Q — )■ this becomes a delta 
function concentrated about Eq.Q. 

3.3.2. Case II 

Here the probability density function is given in terms of the Bessel poly- 
nomials and Whittaker functions (here expressed as 2-fo functions): 

l7i| -^r (-2ri-A-l)^£^„(. 



0<ra<- 
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+ 



27rr(i/x)r(— z/i) 



(26) 



where z = 7ix — 72 and z' = 7ix' — 72 
and 



2-^0 



0-A72 
71 



(27) 



(28) 



Unfortunately, no corresponding simplification of this via a Mehler formula 
is possible. The stationary part that leads to Eq. fl2^ is, 



P. 



{II) 



St 



( \ - l7i| (-/3i-l) _/!2 

^""^ Pt^' r (-A) " ' ^ • 



(29) 



5.5.5. Case /// 

Our efforts for this case have not satisfactorily led to formulae analogous 
to Eq. fl26l) . However the stationary density is straightforwardly derived. We 

redefine z, j3i and (32 for this case as 2; = (tJ"^ ~ {a-nd similarly for 

2;'), and 



Pi 



2 



W7i 



^il ^ ' r]7?v^T^ V 72 

The stationary density is then. 



a 



Ac 



(30) 
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l^^l ^^2 _j_ ]^'j^ig^2 arctan{2) 



(31) 



(cos6»)2(/5i+i) 



These are the solutions of the Fokker-Planck equations, including the 
time-dependence of the probability distribution functions for two out of three 
of the cases. Armed with these, we examine the stability question. 
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4. Stability considerations 

Recall that deterministic stability of the phase synchronised fixed point 
is characterised by the condition that, when initial fiuctuations are small and 
the coupling sufficiently large then under time evolution the system remains 
in the vicinity of the fixed point, Eq.(|9]). Here, if for some modes r this 
condition is not satisfied then deterministic trajectories will drive the system 
away from phase synchronicity. Now, in the presence of noise the system is 
no longer described by a single point but by a diffuse region governed by the 
probability density P{x, t) which is a solution to the Fokker-Planck equation 
Eq.f lT5|) . Using the graph Laplacian, we represent this as a region in the 
mult i- dimensional space with axes Xr through the product form P{x,t) = 
Ylr Prixr,^)- Followiug Schuss [loj , we determine now whether the system 
can diffuse beyond the basin of attraction for some particular mode r, namely 
whether there is a finite probability - as determined from the solutions now 
to Eq.f lT7|l - that 1x^(^)1 > 1. When such modes occur, phase synchronicity 
can be said to be stochastically unstable. Moreover, where we have access to 
the complete time-dependence (and not just stationary distributions) we can 
determine whether some solutions to the truncated equations leave the region 
^ 1 in finite time; in other words, there may be cases of instability that 
the stationary distributions do not detect. 

4.1. Case I 

For the case of 71 = 0, as we can see in the time-dependence of Eq.f E^ . 
the stochastic system tracks the deterministic system which itself is stable, 
never departing from the basin of attraction of the fixed point x* defined in 
Eq.(|8]). We can also see from the stationary part of the probability distribu- 
tion Eq.f l24p that the peak is determined from 

dx - U ^ ^peak - 

We can go further and track the movement of the density peak of the diffusing 
system with respect to time using, 

^pi'){x,x%t') = Q^x = x'e-^^^'-'^ + ^ fl - e-'^^^*-*')) . 
ax a\ \ J 

For small diffusion constant VL the phase synchronicity is therefore effectively 
stable. The tail of the distribution, however, extends outside the region 
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\xr\ ^ 1. But from Eg. (12^ . the probability of the diffusing system escap- 
ing the fixed point is exponentially suppressed. Thus stochastic fluctuations 
on the frequencies alone (71 = 0) do not significantly destabilise phase syn- 
chronicity. 

4.2. Case II 

In this case, we cannot so easily read off the final distribution Eg. (126 p 
whether there is a probability that the system will leave the regime 1x^1 <^ 1. 
In the section below we will plot various cases which will give more insight 
into the evolution of the probability distributions. Analytically, we can again 
locate the peak of the stationary probability distribution from Eq. (!29|) : 

^P^'\x) = ^ a: Hi = ^^ + 7r^ = ^jf\ (32) 
dx '^''''^ 1 -u M 1 -I- M 

We thus see that the peak is shifted compared to the deterministic fixed point. 
However, noting that the denominator of Eq.f l32p is always greater than one, 
the only means for ^^j^^^j^ to exceed one (for a strong instability) is if 7172 

are sufficiently large. In particular, this means that stochastic fluctuations 
of the coupling alone (72 = 0) give 

(72=0) ^ ^ ^ 1 

•^peak 1 I EZl 

In this case, the tail can extend beyond |a;| = 1 but, for sufficiently small 
Q, the probability that the system can drift from the basin of attraction is 
suppressed by a combination of exponential and power-law factors, as can be 
read off the large z properties of Eg. (1290 . 

Evidently, given that increasing 71 will increase the denominator of Eq. fl32l) 
more than the numerator, the parameter 72 is the most effective means for 
pushing the peak of the distribution past x = 1. The coefficient of the 
stochastic fluctuations of the frequencies, while also maintaining fluctuations 
of the coupling, is the effective means of driving the system into instability . 

4.3. Case III 

We can again locate the peak of the stationary probability distribution 
from Eq.dn]): 

(33) 

o-A 
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The noise correlation parameter c becomes an additional degree of freedom 
for driving instability, with c > 1 for fixed 71 and 72. However, even with 
the peak of the distribution lying in the basin x < 1, the tail of the distribu- 
tion is different from cases I and 11. Noting that the leading behaviour in a 
Taylor expansion of e'^^arctanz ^ constant, the behaviour of the stationary 
distribution for large x is purely power-law. 

This power-law dependence for the correlated case means that the suppres- 
sion is considerably weaker than for the previous cases. 

We conclude therefore that increased correlation between noise in cou- 
pling and frequencies is the main driver of stochastic instability from phase 
synchronisat ion . 

5. Probability Density Plots 

When the noise is correlated we have been unable to infer directly from 
the analytic form of the probability distributions whether the system can be 
driven from the basin of the fixed point in finite time, and have thus relied on 
stationary limits. For the cases where we have solved the time- dependence 
we plot several cases to check this possibility. 

To place the stochastic dynamics back into the network context we use 
a case considered by one of us previously in j3] for which synchronisation 
properties of the deterministic system are known. We use a small-scale N = 
27 node variation on the Ravasz-Barabasi (RB) graph [19] that combines 
elements of hierarchy and randomness. In Fig{T] we show the graph in a 
circular embedding and a plot of the spectrum of its Laplacian eigenvalues. 
The four low-lying eigenvalues are Ai = 0.586, A2 = A3 = 0.764, A4 = 0.944. 

We also used, in [3], a set of intrinsic frequencies coi drawn from a uniform 
distribution in the range Ui G [0, 1]; the choice of maximum value of order 
unity fixes a free scale of the model. Figj2] shows the particular frequencies 
used in this case. For these frequencies, the largest values of the ratio pr = 
cu^^yXr (which determines the fixed point) are = 0.530, pi = 0.339, ps = 
0.241. Note that both r = 1,4 correspond to eigenvectors with all non- 
zero components while, for example, u^^^ = (O5, —1.236, — 14, O5, 1.236, 14, O7) 
where x„ here means that the n sequential components have value x. We 
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therefore examine this mode, r = 3, because it involves fewer participat- 
ing nodes. The corresponding projection of the frequency vector onto this 
eigenvector is u^^'' = 0.1843. 

According to first order considerations of the deterministic system 0], 
for all couplings a > p4 = 0.53 the system should be stably locked in a 
state of global phase synchrony. Numerical calculations confirm this 0]. In 
fact, this behaviour extends down to coupling of o" = 0.212 whereupon the 
system exhibits semi-stable periodic behaviour (consistent with the second 
order considerations of 

We therefore know that the deterministic system has a stable phase syn- 
chronised fixed point at a coupling, say, of a = 0.8. We examine the effect 
of noise on the r = 3 Laplacian mode. On this basis, we expect a valid 
fixed point at x* = ps/cr = 0.3. The initial configuration is also taken in 
the vicinity of the fixed point, x' = —0.2 (we use a negative initial value to 
give scope for the probability distribution to detectably shift in the narrow 
interval x G [—1, 1]). We take the diffusion constant to be = 0.1 and use 
an initial delta function at t' = 0. 

We consider three cases, one with purely additive noise and two combining 
aditive and multiplicative noise: 

• Case I: 71 = 0, 72 = —2, for which the Sturm-Liouville operator Eq. (l^ 
has a purely discrete spectrum with the Hermite polynomials as eigen- 
functions; 

• Case II. 1: 71 = 2, 72 = —2, which involves a mixed spectrum of Eq. fl2ip 
and Bessel polynomials; and 

• Case II. 2: 71 = 2,72 = —7, with a similar analytic structure to the 
previous choice. 

Because we do not have access, at this stage, to the time-dependence for Case 
III we shall instead compare the stationary distributions across these cases 
with the given graph structure and frequency spectrum. We thus compare 
the stationary limits of the above three choices with the following correlated 
choices: 

• Case III.l: 71 = 2,72 = —2 and c = 1.1; 

• Case III. 2: 71 = 2,72 = —2 and c = 0.9; and 
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• Case III.3: 71 = 2,72 = -2 and c = 0.2. 



Using p(-f) given in Eq. fl25|) with the following parameter values 
u = 0.184, a = 0.8, A = 0.764, 71 = 0, 72 = -2, 1] = 0.1, a;' = -0.2. 

The t = plot in Figj3]is close to a delta function at x = —0.2, reflecting 
the choice of initial condition. As time progresses the peak of the distribution 
moves to the right as the distribution flattens out, finally reaching its sta- 
tionary shape concentrated around x = ^ = 0.3, reflecting the t — )■ 00 limit 
of Eq.(I7]). For this choice of parameters, the tail evidently extends across the 
boundary x = 1 but the strong (exponential) suppression is also manifest. 
Naturally, smaller values of the diffusion constant Q suppress this further. 

Using P(^^) given in Eg. (12 6 p with the following parameter values 

u = 0.184, a = 0.8, A = 0.764, 71 = 2, 72 = -2, = 0.1, x' = -0.2. 

leads to the evolving probability density shown in FigJH Near t = the 
curve is visibly close to the initial delta function shape around x = —0.2. As 
time progresses the peak of the distribution broadens, becomes asymmetric 
but does not shift significantly. This is because the distribution reaches its 
stationary (t = 00) shape with a peak at x = —0.21, close to its starting point 
and consistent with Eg. (15^ (and Eg. (1551) with c = 1). The distribution is 
strongly suppressed for negative values, consistent with its asymmetry. The 
suppression is less severe for positive values of x, though not easily seen in 
the plot, with still with no significant support for x > 1. 
Using again Eg. fl26|) now with the parameter values 

u = 0.184, a = 0.8, A = 0.764, 71 = 2, 72 = -7, Q = 0.1,x' = -0.2. 

we obtain the evolving density in FigJS] We now observe a distinct shift to 
the left from the initial delta function at x = —0.2, with the density becoming 
asymmetric and settling into the stationary configuration with peak at x = 
— 1.2, consistent with Eg. (132 p . The peak is evidently outside the region 
|x| < 1. Moreover, the tail of the distribution is 'fat', evidently 'spilling' 
across the boundary x > 1 with weaker suppression than before. 

We note overall that the shifts in the distributions as time evolves do 
not show multiple 'bounces' across the interval |x| < 1, rather a monotonic 
evolution from the initial configuration to the stationary distribution. In- 
deed, even with small noise the distributions become too diffuse to allow any 
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discernment of changes of direction of the system within the narrow bounds 
of the region |x| < 1 from subsets of trajectories. Therefore, if the stationary 
density has a peak lying inside the interval |x| < 1 then for all intermedi- 
ate times the support of the distribution will also lie in this domain. Then 
stochastic instability can be concluded from the stationary density. 

Finally, we compare in FigE] the (appropriately normalised) stationary 
densities across all the cases, including those with correlated noise. We 
observe, consistent with the above time- dependent plots and our analytic 
results, that as we move to more and more complex combinations of noise, 
both additive and multiplicative, the position of the peaks of the long time 
distributions move outside the basin of attraction |x| < 1. However, the 
height of the peaks reduces and the degree of smearing increases. We also 
observe a smooth transition between the curves from Case II to Case III. 
Though we do not plot it, increasing c > 1 further suppresses the height of 
the peak and smears the support of the distribution across broad regions of 
values of x. 

6. Summary and Discussion 

We have demonstrated the circumstances under which instabilities can 
be generated by stochastic fluctuations in the Kuramoto model of oscillators 
coupled on a general network. This was achieved by applying both additive 
and multiplicative noise to the dynamical equations linearised in the vicinity 
of the phase synchronised fixed point where the equations can be decoupled 
into eigenmodes of the graph Laplacian of the underlying network. With the 
noise the dynamics is no longer strictly deterministic but is described prob- 
abilistically. In a number of cases, we exactly solved for the time-dependent 
probability density function for each Laplacian mode of the system. How- 
ever, the stationary densities provided sufficient information to determine 
the stability of the system. 

For sufficiently large diffusion constant the tails of the probability distri- 
butions always extend beyond the deterministically defined fixed point basin 
of attraction and provide for a finite probability that the mode can drift from 
conditions consistent with synchronicity. However even for small values of the 
diffusion constant correlation between stochastic fluctuations of the frequen- 
cies and coupling generate a signiflcant mechanism for instability. Firstly, 
suppression of the probabihty outside the basin of attraction is weaker as the 
distribution tail changes from being exponential to power-law. We may call 
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this a weak stochastic instabihty. More pronouncedly, correlation of noise 
allows the peak of the stationary distribution to be shifted outside the basin 
of attraction. In this case, trajectories will escape the basin of attraction 
with probability of order one rendering the system strongly stochastically 
unstable. 

We have thus far avoided limiting our calculations to particular networks 
popularly used in the literature or particular frequency distributions. We 
can now make some broad statements about such cases while treating the 
frequency distribution as fixed. As mentioned at the outset, for deterministic 
systems on easily disconnected graphs larger coupling is needed to stabilise 
phase synchronisation due to the ratio With noise now, we note that in 
determining the peak of the stationary density for generally correlated noise 
with Eq. fl33|) . the stochastic parameter combination cf27| 72^'* is played off 
against the deterministic system combination a\r. Thus larger values of the 
latter combination (strong coupling or high connectivity of the graph) can 
be overcome by larger values of the stochastic parameters to drive instability. 
Conversely, small values of the stochastic parameters suffice to de-stabilise 
easily disconnected graphs. However there is no 'free lunch'. Since low-lying 
modes r ~ 1 for poorly connected graphs (A^ < 1) typically have a larger 
number of nodes in the corresponding eigenvector the noise must be 
manifested in a large number of nodes for the system to be stochastically 
unstable. Depending on the nature of the system and whether stability is 
desired or undesired, it may be expensive to 'attack' or 'defend' many nodes. 

The balancing point of stochastic instability is somewhere in the middle: 
graphs whose distribution of Laplacian eigenvalues exhibits a 'bulk' around 
Af. ~ 1 for 1 ^ r ^ — 1, namely graphs with a large number of subgraphs 
that can be disconnected by a small number of link removals. This is the case 
for sparse Erdos-Renyi and Scale- Free graphs |2o|. Such graphs are therefore 



easily susceptible to stochastic instability. Contrastingly j20[, dense Erdos- 
Renyi (if generated with high rewiring probability from a regular graph) 
and Small World graphs have Laplacian spectra with the bulk of eigenvalues 
shifted away from A^ = 1. These graphs are thus less susceptible to stochastic 
instability. For the Small- World case this applies even for low re-wiring 
probability, because even then the bulk of the Laplacian spectrum is shifted 
away from A^ ~ 1. This mirrors the strong synchronisability properties of 



the Small- World graph observed by [21|. Our analytic results, such as the 



positions of peaks of probability distributions, are sufficiently straightforward 
that the 'weak points' - those nodes significant in determining the transition 
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to stochastic instability - can be easily computed for any network of interest. 

It is difficult to find a precedent for our results in the literature given 
the limited cases that have been investigated. For identical frequencies and 
a complete graph, Park and Kim showed that noisy coupling - or mul- 
tiplicative noise - can induced non-trivial behaviours: a transition of the 
synchronised properties of the system from there being one to two clusters. 
Combining additive and multiplicative noise in the completely connected 
identical oscillator case leads to even richer behaviours, with the former able 
to suppress the effects of the latter as well as inducing saddle-node bifurca- 
tions Q. In our case, we have been able to go further than js], lof in analysing 
stability of any network that can be described by an adjacency matrix. On 
the other hand, apart from identifying conditions for instability we cannot 
say more than these authors about alternative stable or meta-stable states 
(for example we cannot see in Eq. fl33|) a mechanism for mutual cancella- 
tion between noisy coupling and frequencies, as observed by joj). We also 
note that both use additional 'pinning' terms not found in the original 
Kuramoto model. 

To make more specific statements about transitions from phase syn- 
chronicity to alternative fixed points or meta-stable states requires stochas- 
tically probing beyond the first order approximations treated here. For ex- 
ample, in our earlier work [3| we observed that second order approximations 
promoted the 'population model' analogue to the Kuramoto model from the 
logistic to the Lotka-Volterra equations. This offers the tantalising possibil- 
ity of exploiting work on the impact of noise on multi-species competition 
dynamics to make more precise conclusions about phase synchronisation. 
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Appendix A. Properties of the Graph Laplacian 



We repeat here the definition of the graph Laplacian [11] for a graph G: 
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where Aij is the graph adjacency matrix and Dij is the diagonal matrix of 
degrees di = Y^jAij. 

The spectrum of Laplacian eigenvalues A,., r = 0, ... — 1, is positive 
semi-definite 



Ao = < Ai < As... < Xn-i 



as can be proven using the Rayleigh-Ritz ratio for symmetric matrices [11 . 

The degeneracy of the zero eigenvalue corresponds to the number of dis- 
connected components into which the graph can be separated. For a con- 
nected graph forming a single component the eigenvector corresponding to 
the (now) non-degenerate zero mode is easily seen to be 

.<».^^(i,...,i). 



The first non-zero eigenvalue is known as the Fiedler [12[ and encodes im- 
portant properties of the graph. It is also known as the algebraic connectivity 
|ll| and provides a lower bound on the variety of measures of connectivity 
of G: 

Ai(G) < k{G) < A(G) < 6{G) < (di, d) < A(G) (A.l) 

where k(G') is the vertex connectivity (the minimum number of vertices that 
need to be removed to disconnect the graph), \{G) the edge connectivity 
(minimum number of edges to be removed to disconnect the graph), 5{G) 
the minimal degree and A(G') the maximal degree of the graph. 

Mohar 22| summarises a range of theorems about the spectrum of the 
Laplacian and relationships between it and invariant distance measures. Chung 



Lapi 

3, 



23|, using a variation known as the 'analytic Laplacian', proves that the dis- 
tance between any two subsets in G is bounded from above by a function of 
the lowest and highest eigenvalues. 

Appendix B. Deterministic limits of the stochastic system 

Recall that the logistic-like Eq.(|6]) with initial condition, Xr{t') — can 
be solved exactly, with solution 

(r) 

^i-a\r{t-t') , / 1 _ -<TA.{t-t') 



Xr{t) = x;e-'^^'-(*-*^ + V 1 - e-'^^'-^*-*^ . (B.l) 
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We refer to this as the deterministic solution to the stochastic system. This 
result can be arrived at from within the stochastic system using both the 
Langevin and Fokker-Planck equations. 

The exact solution to the corresponding Langevin equation is given by, 



ti=t 



and the integration constant, kr, is given by. 

We observe that the zero limit of all 7 parameters gives back the deterministic 
result Eq.dZ]). 

We can also see the deterministic result arising from the solutions to the 
Fokker-Planck equation as Q approaches zero. This statement can be ele- 
gantly demonstrated in the pure stochastic frequency case, where the prob- 
ability density is given initially through the Hermite polynomials. Starting 
with the expression Eq.f l25p for the density after use of the Mehler formula 
and applying the following delta function identities, 

1 f x'^\ 1 
6{x) = lim exp I , 6{ax) = —6{x) , a G M 



we obtain, 

lim P (x, x'lt, i!\ = b(x- x'e"'"^* - — f 1 - e~'"^* 



The result is a pure Dirac delta function about the deterministic trajectory. 
We also note that 



(x) = x'e-"^* + — ( 1 - e-"^* 
crA V 

so that the expectation value of x tracks the deterministic trajectory. 
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Appendix C. Spectral categories for Fokker-PIanck equations 



Following the work of Linetsky 2J], we detail the form of g{x, t) which are 
solutions to Eq. fl2U]) . the remaining part of the Fokker-PIanck equation after 
separation of the stationary solution. This solution depends on the detailed 
form of the functions q{x) and s{x) and the range of x. 

Consider the Sturm-Liouville operator Eq. (!2T!) . and the associated eigen- 
value problem Eq. (!22!) . namely 

s{x)—u{x) + q{x)—u{x) + ^u{x) = , ^ > (C.l) 

where the range of x for Eq JC.ll is given hj x E [xmin,Xmax]- Let e be an 
end-point of this range, namely e G {xmin,Xmax}- 

If e 7^ ±oo then Eq. llC.ll) is classed as non-oscillatory for all ^. If e = 
±oo, refered to as a natural boundary, then we must examine more closely. 
Transforming the variable x and the solution u{x) via the following, 

yix) = [ , viy) = u{xiy))[Wixiy))]^sixiy)t^ (C.2) 

the corresponding Sturm-Liuoville equation for v{x{y)) becomes the one- 
dimensional Schrodinger equation with the potential. 



^/six) 

The classification at the natural boundary is then given by the following 
theorem. 

Theorem Classification at natural boundaries (Linetsky) . 

• If e is transformed into a finite endpoint under Eq. (1C.2l) . namely y{e) ^ 
±oo, then Eq. llC.lll is classed as non-oscillatory for all C,. 

• If the transformed end-point is infinite, y{e) = ±oo, and lima,_>e V{x) = 
{oo, 0}, then Eq. OC.ip is classed as non-oscillatory for all ^. 

• If for some finite non- vanishing V{e), we have lim^.^e^(3;) = V{e), 
and \imx^ey'^ix){V{x) — V{e)) > — ^, then Eq. flC.ip is classed as non- 
oscillatory for ^ G [0, V{e)] and oscillatory for ^ > V{e). 
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We are now equipped to classify the exact forms of the solution to Eq. ( 12U]) . 
More detail on these aspects can be found in section 4 of 24 . 



Appendix C.l. Spectral category I 

If both boundaries exhibit no oscillatory behaviour, then the spectrum is 
purely discrete and the probability distribution is given by, 

oo 

P{x, x'\t, t') = Pstix) J2 e"""^*-*' Vn<^n(x)<^.(a;') , p„ G M (C.4) 

n=0 

where (pn{x), n > 0, are eigenf unctions according to 

d 

^i^)l-^'^n{x) + q{x)—ipn{.x) + ?7„(/)„(x) = (C.5) 
ax^ dx 



for the discrete eigenvalue ?7„. Following Chap. 5 of [15[, the spectrum is 
positive semi-definite: 

?7„ e M , = r/o < r/i < < • • • • 

The eigenf unctions are the classical orthogonal polynomials, for specific 
forms of s{x) and ranges of x: 

• If s{x) is a constant and x G M then V9„(x) takes the form of a Hermite 
polynomial. 

• If s{x) is linear and x G M+ then (pn{x) takes the form of a Laguerre 
polynomial. 

• If s{x) is quadratic with two distinct real roots and x G [—1, 1] then 
<y9„(x) takes the form of a Jacobi polynomial. 

These classical orthogonal polynomials come with the following infinite or- 
thogonality property, 



dxPstix)(pnix)(pmix) = dmnK < OO for m, n = {0, 1,2,.. .} (C.6) 

n 

indicating that the infinite set of eigenfunctions are square integrable with 
respect to the weight function Pgt- In Appendix F we use the above or- 
thogonality, and the initial condition P{x,x'\t' ,t') = 6{x — x') to obtain 
normalisation values, p„, in Eq. flC.4p . 
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Appendix C.2. Spectral category II 

If one of the boundaries exhibits no oscillatory behaviour, with the other 
boundary exhibiting oscillatory behaviour for C, > V{e), then the probability 
density is given by, 

N 

P{x,x'\t,t') = Pst(x) [^e-''"(*-*'Vn</^n(a;)y^„(x') 

n=0 

POO 

/ ci/ie-^('^)(*-*'V(/x)V'(/i,a;)V'(/^,a:') • (C.7) 







We treat the discrete and continuous parts of this separately in the following. 

Appendix C.2.1. On the discrete eigenvalues 

For the finite sum term in Eq. (lC.7|) . the eigenfunctions '^n{x) satisfy a 
differential equation equivalent to Eq. (lC.5p . with a corresponding discrete 
eigenvalue ?7„, but where the orthogonal polynomials no longer exhibit the 
infinite orthogonality property Eq. (IC.6p . but a finite orthogonality property, 
namely 

nax 

dxPst{x)iPn{x)frn{x) = 6,nnhn < OO m, U = {0 , 1, 2 , . . . , N} (C.8) 

.71 

for finite integer A^. We refer to these as non classical orthogonal polynomi- 
als. The cases of s{x) which generate these polynomials are as follows: 

• If s{x) is quadratic with one repeated root and x G M+ then ipn{x) 
takes the form of a Bessel polynomial. 

• If s{x) is quadratic with two distinct real roots and x G M+ then V9„(x) 
takes the form of a Fisher-Snedecor polynomial. 

Appendix C.2.2. On the continuous eigenvalue 

For the continuous part of Eq.f lC.7p we have the following eigenvalue 
equation for ip{^,x), 

d'^ d 

x) + q{x)—tlj{n, x) + A(/i)V'(/i, x) = (C.9) 

for the continuous eigenvalue A(/i). The integration variable /i provides a 
continuous index for the eigenvalue A(yu) (as n provides for the discrete value 
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rin). Again, following Chap. 5 of 15| we have the following condition for the 
eigenvalue, 

A{fi) e M , A(yu) > 

The eigenfunction, ilj{fi,x), is a fundamental solution to Eq. (IC.9P which is 
valid in the neighbourhood of the non-oscillatory boundary. The eigenfunc- 
tions for continuous eigenvalues must be orthogonal to eigenfunctions for the 
discrete rjn, n E {0, . . . , N}. 

The eigenvalue A(/i) can be found in one of two ways, from the finite 
orthogonality property or simply from considering the oscillatory cutoff point. 
Using the finite orthogonality approach (for the other method see [2J]) we 
use the property Eq. (]C.8l) . to infer the existence of a real number, J', which 
is greater than the largest discrete eigenvalue, r]^, 

Vn<J. (C.IO) 

If this is not the case then the finite orthogonality property (Eq JC.Sp fails, 
leading to the catastrophic condition (specific details of this statement will be 



given in [Appendix Fj), /;,„ — )■ oo. Thus we define the continuous eigenvalue 



in the following way, 

Aifi) = J + Kfi\ {k, fi} > 0. (C.ll) 

We have necessarily that J' = V{e). The choice of /x^ for the dependence 
on the index /i in this parametrisation guarantees that the eigenvalue ranges 
from [V^(e), oo]. 

The eigenfunctions with continuous eigenspectra are generally not poly- 
nomials but obey an equivalent orthogonality condition amongst themselves, 



dxPst{x)ip{fi, x)ip{i^, x) = 6{fi — i^)h{fi) (C.12) 
and are completely orthogonal with their discrete spectrum counterparts, 

dxP,t{x)^{^i,x)'fn{x) = 0,n={0,l,2,...,N} (C.13) 



As with the first spectral categorty, the orthogonality conditions and the 
initial condition P{x,x'\t',t') = 6{x — x') are used to obtain the normal- 
isation constants {p„,p(/i)}. Specific instances of these will be solved in 
Appendix F 
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Appendix C.3. Spectral category III 

If both boundaries exhibit oscillatory behaviour for some ^ > oo), V{oo)} 
and V{—oo) = V{oo) then the probability density is given by, 



Pix, x'\t, t') = P,tix) e"''"^*~*'Vn<^n(a:)v^„(x' 

n=0 



+ 



2 

E 

ij=i 



(C.14) 



Note that there exist more general versions of spectral category III where 
V^(— oo) 7^ V{oo), but we do not need to consider these in this work. 

Appendix C.3.1. On the discrete eigenvalues 

As with the previous spectral category, the eigenfunctions V5n(x) in Eq. (lC.14p 
for the discrete eigenvalue ?7„ satisfy an equivalent differential equation to 
Eq. flC.SP and exhibit a finite orthogonality property Eq. flC.Sp . The relevant 
restriction for s{x) for this spectral category is, 

• If s{x) is quadratic with two distinct complex roots and a; G M then 
(y9„(x) takes the form of a Romanovski polynomial. 

The continuous part of the spectrum continues to present some diffulties 
which is why we do not consider this case any further in this paper. 



Appendix D. Properties of special functions for solving Fokker- 
Planck equations 

In this section we give basic properties of the hypergeometric functions 
and differential equations used in our calculations. To begin we give some 
necessary definitions. 

Appendix D.l. The Pochhammer symbol and Gamma functions 

Two quantites we call upon extensively in the definition of hypergeometric 
functions are the Pochhammer symbol, (x)„, and the Gamma function, T{x), 
given by, 

ix)n = + •^■) = ' r(^) = / dte'H^-' , M(x) > 0. 
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Appendix D.2. Hypergeometric differential equations 

In this work we consider hypergeometric differential equations of the form, 

d'^ d 
s(x)—^ilj(fi, x) + q(x)—i'(fi, x) + A(fi)ip(fi, x) = (D-1) 
dx^ dx 

where, 

s{x) is at most a quadratic polynomial of x 

q{x) is a linear polynomial of x 

A(/i) the eigenvalue - is not a function of x. 

Appendix D. 3. Hypergeometric functions 

Solutions to the above differential equation are referred to as hypergeo- 
metric functions, generally given as. 



p-Tq 



CtX ... dp 
hi ... bq 



j=0 



In this work we are limited to p = {1,2} and q = {0,1}. For more information 
see Chaps. 13, 14 and 15 of [l6|, and references therein. 

Appendix D.^. Order n polynomials 

Through the following choice for the eigenvalue, 

f d 1 d'^ \ 

Aifi) =Vn = -n {-^q{x) + ^{n - l)—s{x)\ (D.3) 

the value of ai in Eq iD.2l becomes — n, and hence the series terminates at 
j = n, leading to a polynomial solution, fn{x), of order n. 

Appendix D. 5. Weight functions 

Each set of polynomials come equipped with a weight function, given by 
Pearson's differential equation, 

^ {s{x)Prmst{x)} - q{x)Pst{x) = ^ Pstix) = CXp | j ^Xj^^ . 
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Appendix D. 6. Canonical forms 

The following lists detail the canonical forms of s{x), q{x) and the weight 
functions of the six distinct cases of the previous section, 

Polynomial s(x) q(x) Eigenvalue 

Her mite 1 —2x 2n 

Laguerre x P + 1 — x n 

Jacobi 1 - Pi- P2- + P2 + 2)a; n{n + /3i + /^s + 1) 

Bessel x"^ + 2)x + P2 -n{n + /3i + 1) 

Fisher- Snedecor x{x + 1) {/3i + (32 + 2)x + /3i + 1 -n{n + /3i + + I) 

Romanovski x^ + 1 2(/3i + l)x + [32 -n{n + 2/3i + 1) 



Label 




Interval 


Pst(x) 




Restrictions 






X G M 












X G M+ 


x^e~^ 




/3 > 


p^/9i,/32)| 




X G [-1,1] 


(l+x)''^( 


1 - xf^ 


{/3i,/32}>-l 






X G M+ 


x'^^e 




/3i <0,/32 >0,n(^) 






X G M+ 


x^^il + x; 




/3i > -l,n(^) 






X G M 


(x2 + 1)^^' 


g/32 arctan x 


/3i < 0,n(-f^) 



Appendix D.6.1. On the restrictions 

The interval of x refers to the region of integration necessary to satisfy 
the orthogonality conditions. Additionally, the restrictions n^^\ n^^'^ and 
for the Bessel, Fisher-Snedecor and Romanovski cases respectively are, 

{n(^),n(^),n(«)} ^ |n < < < - (/^i + ^ 

(D.4) 

These two conditions allude to the finite orthogonality conditions held by 
these sets of polynomials. We shall detail the derivation of these conditions 
shortly. 

Appendix D. 7. Rodrigues formula 

Aside from solving the differential equation Eq. (IC.SP directly, it is possible 
to obtain the explicit form of each order n polynomial using the so called 

Rodrigues formula, 

1 (i" 

^n{x) = C„— --— — {P3t(x)s"(x)} where C„ G M. (D.5) 

l^st\x)s[x) ax" 
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Appendix D. 8. Orthogonality 

We now detail the orthogonality of the polynomials. To do this we shall 
explicitly derive the infinite Eq. flC.6p . and finite Eq.f lC.Sp orthogonality of 
the classical and non-classical polynomials respectively. 

We begin by recasting Eq. llC.SP into the following form, 



d 

dx 



s{x)Pst{x)^(fn{x] 



+ T]nPst{x)ipn{x) = 0. 



(D.6) 



Consider two copies of the above equation, one with index n and one with 
m, m n. Multiply equation n with iprn{x), and multiply equation m with 
(Pn{x). Finally, we subtract the two corresponding equations to obtain. 



Pstix) ipnix) Lprnix) 



1 



Vn{x)-^ \s{x)Pst{x)-^iPm{x) 

dx dx 



d 

dx 



s{x)Pstix)-^(Pnix) 



Integrating both sides with respect to x over the relevant intervals, and per- 
forming integration by parts to the right hand side we obtain. 



dxPst{x)(pn{x)iprn{x) 



1 



d 



s{x)Pstix) <j (pn{x) [ -^^rnix) 



- 1 ) ^m[X) 



(D.7) 



Appendix D.8.1. Infinite orthogonality 

To illustrate the infinite orthogonality property of the classical orthogonal 
polynomials, consider the example of the Hermite polynomials. Inserting 
'^ni.x) = Hn{x) into the right hand side of Eq. flD.7p we obtain. 



le \ Hn{x) ( ^Hm(x)] - ( ^HJx) 1 H^x) 



dx 



dx 



order m + n — 1 polynomial 
We notice that the order m+n — 1 polynomial is completely suppressed by the 

2 

weight function, , in the x — ?■ ±oo limit for all m and n, hence showing 
the infinite orthogonality of the Hermite polynomials. It is an equivalent 
process to show that the remaining classical polynomials have this property. 
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Appendix D.8.2. Finite orthogonality 

To illustrate the finite orthogonality property of the non-classical orthog- 
onal polynomials, consider the example of the Bessel polynomials. Inserting 
(x) = Bn^'^^\x) into the right hand side of Eq. flD.7l) we obtain, 



\x^0 ■ 



order /3i + m + n + 1 polynomial 

Recall that /5i < 0. Again, the polynomial part in the above expression is 
of order /3i + m + n + 1, and hence is completely suppressed by the function 
e ^ , in the a; — )■ limit for all m and n. However in the a; — )■ oo limit, e ^ 
does not suppress the polynomial expression. Hence in order for the limit 
to be zero, we require m + n < —{/3i + 1). Performing an equivalent power 
counting argument for m = n, in order to satisfy the following bound on the 
orthogonality integral. 



POO 

/ dxe-^x^^{Bi[^^'^'\x)) <oo 
Jo 



2 



we obtain the already given restriction of n, n < 

It is an equivalent process to derive the corresponding restrictions as- 
sociated with the Fisher-Snedecor and Romanovski polynomials. For more 
information on the classical orthogonal polynomials (Hermite, Laguerre and 
Jacobi), see Chap. 22 of [16]. For more information on the Bessel poly- 



nomials see 17, 25 and references therein. For more information on the 



Fisher-Snedecor polynomials see j26l| . Finally, for more information on the 



Romanovski polynomials see |l8l . l27l . |28 . 



Appendix D. 9. A word on the Fisher-Snedecor polynomials 

The Fisher-Snedecor polynomials are actually the Jacobi polynomials un- 
der a scale/translation transformation and greatly increased range. It is this 
increased range which causes the finite orthogonality property, hence forcing 
us to place the Fisher-Snedecor polynomials in the non-classical category. 



Appendix E. Inapplicable parameter cases 

Appendix E.l. Laguerre polynomials 

For this case we require that s{x) be linear. However, careful inspection 
of s{x) in Eq. ( 1T6|) reveals that there is no sensible method to obtain a linear 
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s{x). Although not apphcable in this context, those interested in applying 
Laguerre polynomials to Fokker-Planck type equations in other contexts are 



referred to [29| and references therein. One obtains a result similar to the 
Hermite case, where the infinite bilinear sum of Laguerre polynomials is 
drastically simplified to a single modified Bessel function of the first kind 
using the Hille-Hardy formula. 

Appendix E.2. Jacobi and Fisher- Snedecor polynomials 

For this case we require that s{x) be quadratic and contain two distinct 
real roots. However, like the Laguerre case, careful inspection of s{x) reveals 
that there does not appear to be a sensible way to achieve this. One could 
attempt to make certain variables complex, but this would have negative 
consequences on the corresponding Langevin equation. Although not appli- 
cable in this context, those interested in applying Jacobi and Fisher- Snedecor 
polynomials to Fokker-Planck type equations in other contexts are referred 
to 



29[ | and jSOj respectively. 



Appendix F. Detailed derivations of spectral cases 



We now begin to assign values to the parameters, {71, 72, c}, in Eq. flC.SP 



This allows us to generate a specific form for s{x), (constant, linear or 
quadratic), hence allowing us to pin down the form of the eigenvalues and 
eigenfunctions of the various cases. 

Appendix F.l. Parameter values - Hermite polynomials 

For this case we require that s{x) be a constant, hence we assign, 

7i = and {72, c} G R ^ s{x) = (F.l) 
The corresponding Langevin equation is 



X = u - a\x + 72cri + 72 V |1 - c \^2- 

Appendix F.2. The stationary density /weight function 
Substituting the above values into Eq. (fT8l) we obtain, 

^st(a;) = 7^ exp <j I dx[x-^ 

V 





\ 2a\ 




[ ^ll 











x-^] } (F.2) 
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where we have assigned the value (^)^ for the integration constant in order 
to complete the square. 

Appendix F.3. Eigenvalue equation 

Substituting the parameter values in Eq. (IF.ll) into Eq. (10.51) and perform- 
ing the following change in variables, 



X 



we obtain, 



d? d 2 

— LPn{z) - 2Z — Lpn{z) = rr]n^n{z) = -2nLpn{z) 

az'^ az aX 



(F.3) 



(F.4) 



The eigenvalue, rjn = crXn, is obtained from Eq.( 1D.3|) . The above equation 
is known as the nth order Hermite differential equation whose solution is 
the nth order Hermite polynomial, fn{z) = Hn{z), which has the following 
hypergeometric form. 



HoJz) 



-1] 



,(2n)! 



n! 



H2n+l{z) 



-n 
1 

2 

— n 

3 
2 



or equivalently from the Rodrigues formula 

HJz) = 



dz' 

Appendix F.4- Analysis at the natural boundaries 

Using the details from Appendix C we have the following form for the 
potential V{z) in Eq.dU^]), 

V{z) = — (z^-l)^ lim V{z) = oo. 

2 2— >±00 

Hence we determine that the behaviour of Eq. (]F.3P at the boundaries is 
strictly non-oscillatory, and hence the solution to the probabihty density is. 



P{z,z\t,t') = e-^'y^e-'^^<'-''^pMz)HrXz') 

^72 ^ 



(F.5) 



n=0 
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Appendix F.5. Initial conditions and orthogonality 

We now determine the normalisation parameter, through the initial 
condition, 

P{z, z% t') = 6 (^]l^l2{z - z')^ = ^J^^^i^ - 
and applying the Hermite polynomial orthogonality condition, 

/oo 
dze-^"HUz)H^iz) = 5^„2"n!v^. (F.6) 
'OO 

We obtain. 



'^7c^A I72I 

Pn 



Thus the probability density is 



\ 00 



Appendix F.6. Mehler's formula 
Applying Mehler's formula, 



(F.8) 

we can eliminate the summation altogether from Eq. flF.7p . Expressing the 



solution in terms of the original x variable the final form of the probability 
density becomes. 



as presented in the main body. 
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Appendix F. 7. Parameter values - Bessel polynomials 

For this case we require that s{x) be quadratic and contain one real double 
root, hence we have the two possibilities, 

c= 1 and {71,72} G M ^ = f (71X - 72)' (F.IO) 
72 = and {71, c} G M s{x) = ^x^. (F.ll) 

For this section we shall concentrate on the choice Eq.l lF.lOll . for which the 
case c = — 1 can be obtained by reversing the sign of 72. It is an equivalent 
process to consider the second choice. The corresponding Langevin equation 
is 

x = uj- aXx + (72 - 7il/)ri. 

Appendix F.8. The stationary density/weight function 
Performing the following change in variables, 

x = — + ^ (F.12) 

7i 7i 

and substituting Eq. f lF.10|) into Eq.( IT8|) we obtain. 



where. 



(F,13) 

Qjf Q-fi \ 71 y 



Appendix F.9. Eigenvalue equation - discrete spectrum 
Substituting Eqs. flF.10yF.12p into Eq. flOSj) we obtain, 



(p d 2 

z'^—(Pn{z)+{{f3i + 2)z + /^a} -^'^n{z) = --^VnVniz) = n{n+/3i + l)^n{z) . 



(F.15) 

The eigenvalue, r]n = —^^n{n + /3i + 1), is obtained from Eq.( ]D.3l) . The 
above equation is known as the nth order Bessel differential equation whose 
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solution is the nth order Bessel polynomial, ^n{z) = Bn^^'^^\z), which has 
the following hypergeometric form [3l[, 



or equivalently from the Rodrigues formula, 



n \ / on 



P2 rf" 



PI 



e ^ 



dz""- 



Appendix F.9.1. The highest eigenvalue 

We note that the highest eigenvalue, rjjq, is bounded by the following real 
number, 

^^7? + 1)' 



(F.16) 



as given by Eg. f lD.4p . 

Appendix F.IO. Analysis at the natural boundary 



Using the details from section Appendix C, we have the following form 
for the potential V{z) in Eq. flC.3p . 



V{z) 



log,(x) lim y'^{z){V{z) - V{oo)) = 



Hence we determine that the behaviour of Eq. (IF.ISP at the natural boundary 
is non-oscillatory for rjn G 0, + 1)^ and oscillatory for rjn > + 

1)2. Hence the explicit form of the probability density is given by. 



P(z,z'\t,t'^ 



2 _fh 



J2 e^"("+''i+i)(*-*')p„5i'^^''^^)(^)Si''i'^^)(/) 

(F.17) 



n 



0<n< 



+ j dfxe-^^f'^^'-''^pM{fx,z)ij{fi,z')} . 
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Appendix F.ll. Eigenvalue equation - continuous spectrum 
Substituting Eqs. flF.10yF.12p into Eg. ( 109]) we obtain, 



d 



z^— z) + + 2)z + (3^} —^{fi, z) + ^A(/i)V;(/x, z) = 0. (F.18) 



The canonical form of the Bessel equation in |3lj, (which has poles at z 
{0, oo}), is given as 

z^—tjj{u, z) + {2az + 6} -i-^(z/, z) - {u + a){u - a + l)V^(z/, z) =Q 



dz 



dz 



with solution, ifj^u^z) = 2-fo (a — z^ — l,a + z^|— f). Comparing this result 
with Eq.f lF.lSj) we have 2a = (3i + 2, b = (32, and hence. 



/3i 



Knowing that the above eigenvalue is restricted by Eq. flF.16|) . we let u 
ijj, — \, and hence. 



2-^0 



n^l / (A + 1)' 



2 



+ /iM , /i > 



iji, — - — h ijj, 



z 
T2 



2 ^'2 

which is our fundamental solution in the neighbourhood of the finite bound- 
ary. 

Appendix F.11.1. The Whittaker function 

Following Chap. 13 of jl6[, we have the following relation. 



ip{fi,z) 



T{-2ifi) . r(2z/i) 



(F.19) 



where M(±/i, z) 



z 



/3i + l 



1 ±2i/i 



2; 



The above solution is proportional to the Whittaker junction of the second 
kind of imaginary order 32l. l33|. 



- 2 V Z 



A f 132 

e 2z 



ft. 

2 



(F.20) 
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Explicity the Whittaker functions of the first and second kind, M^^ii^(z) and 
M^K,i//(^) respectively, satisfy the Whittaker differential equation, 

^,W.M + + ^ - 1) W.,,iz) = 0. (F.21) 

The purpose of estabhshing this connection between ip{z, ix) and WK,i^j.{z) will 
become clear when we consider the orthogonality conditions of the continuous 
spectrum eigenf unctions. 

Appendix F.12. Initial conditions and orthogonality 

We now determine the nomalisation factors, p„ and p(//), through the 
initial condition, 

P{z,z'\t\t') = S (^^(^ - z')^ = \^i\6{z - z') 

and applying the following discrete and continuous orthogonality conditions, 

r dze-'^z^^Bt''^){z)Bl^^^^^\z) = (F.22) 

Jo 2n + Pi + l 

(F.23) 

r<x> 

/ dze-^z^'Bl^''^'\z)ilj{ix,z) = (F.24) 
Jo 

where < m, n < — for the discrete case(s). Using the above relations, 
we obtain for p„ and 



Pn = TT-rP'. 



f^|7i|^-/3i-i 2n + /ji + l 
2n! r(-n-^i) 



^ l]|7i| r(^ + ./.)r(^-.;.) 

2 27r/3fi+^r(z^)r(-z^) 

Hence we obtain the complete form of the probability density for this case, 
as given in the main body. 



0<n<— ^. 
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where z = 71Z — 72 and 2' = 7ix' — 72. 



(F.25) 



Appendix F.12.1. A note on the orthogonality relations 

The discrete orthogonahty relations Eq. flF.22p . were given in 17|, and 
the authors noted that the proof requires nothing more than the knowledge 
of how to apply integration by parts. Our continuous orthogonality relation 
Eq.f lF.23"|) . was adapted from the following relation, 



°° dx 



_ 27rr(i/i)r(-i/i) 
{,^,(/}>0 (F.26) 



originally given by Wimp [32j, which is a particular case of a more general 
relation involving the Meijer's G-function. The above condition was rederived 
in 331] by the more elementary means of directly using the properties of the 
hypergeometric functions 2-^0 and iFi. We shall adapt the proof in [s^ 
in the next section for the hypergeometric function 2-^1, thus deriving new 
orthogonality relations. 

Appendix F.13. Parameter values - Romanovski polynomials 

For this case we require that s{x) be quadratic and contain two distinct 
complex roots, hence we assign. 



c 7^ ±1 and {71,72} e 



(F.27) 



The corresponding Langevin equation is 

aXx + (72C - 7ix)ri + 72a/|1 - c^lFs 



X 
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Appendix F.14- The stationary density /weight function 
Performing the change in variables 

72 



X = 

71 



(^Vll -c2| +c) (F.28) 



and substituting Eq. (lF.27p into Eq. lITS!) we obtain, 

2 



PJz 



fi7i(|l-c2|)(^2 + i) 



^^7f 7 \ V|l-c2| V 72 / J + 1 

^ -^2 ^ ^^|/3lg/32arctan{^) (^-p 



fi7i(|l-c2|) 



where, 



/3i = -7^ - 1 , /32 = , , I , -aAc) . (F.30) 

^7i (]7iV|l-c2| V 72 / 

Appendix F.15. Eigenvalue equation - discrete spectrum 
Substituting Eqs. flF.27yF.28p into Eq. flU3|) we obtain, 

(p d 2 

{Z^ + 1} -^^Pn^z) + {2(/3i + \)Z + /32} J-^Vn{z) = —^Vn^niz) 

= n{n + 2(3i + l)^n{z). (F.31) 

The eigenvalue, rjn = —^^n{n + 2/3i + 1), is obtained from Eq.( 1D.3|) . The 
above equation is known as the nth order Romanovski differential equation 
whose solution is the nth order Romanovski polynomial, (pn{z) = R^^'^^\z). 
It is possible to transform the nth order Romanovski differential equation 
into the Jacobi differential equation, given as, 

{l-e} ^¥^n(0 - {A - /32 + + /32 + 2)0 -^MO 

+n(n + /3i + /32 + l)</^n(0 =0 

where 

MO-Pn iO - ^, 2^1 ^^^^ 

(F.32) 
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Performing the transformation z = in Eg. (lF.3ip we obtain, 



+ 2/3i + l)(^„(0 = 0. (F.33) 



Thus, comparing Eq. flF.33p with Eg. (IF. 32]) we obtain the following hyperge- 
ometric form for Ri^^'^''^\z), 

^,«r (ft -4.1)^^ ;-4^^;;^ |i±i£). 

Equivalently, from the Rodrigues formula we obtain, 

Jn 

JliPi'M(z) = (z"^ + l)-/^lc-fearctan(2:) " // 2 , -|^y«+/3ig/32 arctan(2:)\ 

We have not succeeded in progressing the analysis satisfactorily beyond 
this point. 
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Figure 1: An "RB"-type network represented in a circular embedding shown in the upper 
panel while its Laplacian spectrum is shown in the lower panel with the eigenvalues plotted 
in increasing order, from right to left 
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Figure 2: Choice of intrinsic frequencies used in this numerical study which is drawn from 
a uniform distribution between [0, 1] 




t = 0.015 



Figure 3: Plot of the probability density, the solution to the Fokker-Planck equation Ea. l|17| l. for different 
values of time starting from an initial delta function at x = —0.2, with noise parameters Q = 0.1 and 
7i = 0, 72 = —2 which lead to a Hermite decomposition. The series of curves describe the evolution of 
the r = 3 Laplacian eigenmode of the RB graph with uj = 0.184, cr = 0.8, A = 0.764. 



46 




Figure 4: Plot of the probability density, the solution to the Fokker-Planck equation Eq.JTTJ, for different 
values of time starting from an initial delta function at x = —0.2, with noise parameters f2 = 0.1 and 
c=l,7i = 2,72 = — 2 which lead to a Bessel polynomial decomposition. The series of curves describe the 
evolution of the r = 3 Laplacian eigenmodc of the RB graph with ui = 0.184, a = 0.8, A = 0.764. 




Figure 5: Plot of the probability density, the solution to the Fokker-Planck equation Ea. l|17| l. for different 
values of time starting from an initial delta function at x = —0.2, with noise parameters f2 = 0.1 and 
c=l,7i = 2,72 = — 7 which lead to a Bessel polynomial decomposition. The series of curves describe the 
evolution of the r = 3 Laplacian eigenmode of the RB graph with uj = 0.184, cr = 0.8, A = 0.764. 
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